Case study on migration of Anser Fabalis
2020-11-21
1 Overview
This notebook contains supplementary material for the paper:
Predicting spatio-temporal distributions of migratory populations by combining citizen science data and Gaussian Process modelling (2020). Antti Piironen, Juho Piironen and Toni Laaksonen. Submitted
The codes below show how to repeat the steps for fitting the Gaussian process (GP) model discussed in the paper. Fitting the model using the package gplite is very easy, and most of the code below is actually related to manipulating the data and visualizing the results.
2 Load the data
First load the required R-packages
library(chron)
library(rprojroot)
library(gplite)
library(ggplot2)
library(ggspatial)
library(maps)
ROOT <- rprojroot::find_root('.gitignore')
source(file.path(ROOT, 'scripts/subplot.R'))Define a few auxiliary variables and functions.
# these variables will define time ranges that are used to select data
# that are representative about the autumn and spring migrations
SPRING <- list(
# note: year does not matter here
start = chron('1.3.2000', format='d.m.y'),
end = chron('31.5.2000', format='d.m.y')
)
AUTUMN <- list(
# note: year does not matter here
start = chron('1.8.2000', format='d.m.y'),
end = chron('30.11.2000', format='d.m.y')
)
within_interval <- function(date, start=NULL, end=NULL, ignore_year=TRUE) {
if (ignore_year) {
year <- 2000 # just some year, does not matter which
mdy <- chron::month.day.year(date)
date <- chron(paste(mdy$day, mdy$month, year, sep='.'), format='d.m.y')
mdy <- chron::month.day.year(start)
start <- chron(paste(mdy$day, mdy$month, year, sep='.'), format='d.m.y')
mdy <- chron::month.day.year(end)
end <- chron(paste(mdy$day, mdy$month, year, sep='.'), format='d.m.y')
return(date >= start & date <= end)
}
return(date >= start & date <= end)
}This will load the data. Here variable SEASON will choose which season is being analyzed (it should be either 'autumn' or 'spring').
# load the data
SEASON <- 'autumn' # 'spring' or 'autumn'
filename <- 'anser_fabalis_2011-2019.csv'
dat <- read.csv(file.path(ROOT, 'data/', filename), sep = ',', header=T)
dat$date <- chron(as.character(dat$date), format='d.m.y')
# filter only observations that fit the defined spring or autumn days
rows <- sapply(dat$date, function(date) {
if (SEASON=='spring')
return(within_interval(date, SPRING$start, SPRING$end))
else if (SEASON=='autumn')
return(within_interval(date, AUTUMN$start, AUTUMN$end))
else
stop('Unknown season: ', SEASON)
})
dat <- dat[rows,]
dat <- dat[order(dat$date),]
# target variable
n <- dat$fabalis + dat$rossicus
y <- dat$fabalis
# predictor variables (input)
time <- as.numeric(as.times(dat$date)) # count only days
input <- as.data.frame(cbind(time, dat$x, dat$y))
colnames(input) <- c('time','x1','x2')3 Fit the model
The code below will specify the model structure. Here we use initial guesses for the hyperparameters that are relatively close to the optimal ones, in order to make the optimization faster. In practice, the optimization should converge also with less accurate initialization, but will take longer time. (Notice though, that it is possible that with a very bad initialization the optimization will not converge, or it converges to some suboptimal local minimum).
cf_spatial <- cf_nn(
c('x1','x2'),
magn = 3.7,
sigma0 = 4,
sigma = 3.5,
prior_sigma0 = prior_half_t(1),
prior_sigma = prior_half_t(1)
)
cf_temp1 <- cf_periodic(
'time',
period = 365,
prior_period = prior_fixed(),
cf_base = cf_nn(
magn = 1,
sigma0 = 0.5,
sigma = 5.5,
prior_magn = prior_fixed(),
prior_sigma0 = prior_half_t(1),
prior_sigma = prior_half_t(1)
)
)
cf_temp2 <- cf_sexp('time', lscale=290, magn=1, prior_magn = prior_fixed())
cf_temporal <- cf_temp1 * cf_temp2
cfs <- cf_spatial * cf_temporal
gp <- gp_init(
cfs,
lik = lik_betabinom(phi=5.5),
method = method_fitc(num_inducing=200, bin_along='time', bin_count=9),
approx = approx_laplace()
)Optimize the hyperparameters. If you would like to actually run this part (will take some time), set the flag fit_model = TRUE. The model fitting is disabled by default to make this notebook compile fast. For plotting the results, we load a prefitted model by default. Note: if the program gives a warning that the optimization has not converged, you can simply run this cell again (if will start the optimization from were it ended up previous time).
fit_model <- T
if (fit_model) {
gp <- gp_optim(gp, input, y, trials=n, tol=1e-6)
}## Optimizing parameters
## p1: log cf_periodic.sigma0
## p2: log cf_periodic.sigma
## p3: log cf_sexp.lscale
## p4: log cf_nn.sigma0
## p5: log cf_nn.sigma
## p6: log cf_nn.magn
## p7: log lik_betabinom.phi
##
## p1 p2 p3 p4 p5 p6 p7 Energy Iteration
## -0.69 1.70 5.67 1.39 1.25 1.31 1.70 2975.12 0
## -0.13 1.70 5.67 1.39 1.25 1.31 1.70 2975.61 1
## -0.69 2.27 5.67 1.39 1.25 1.31 1.70 2977.06 2
## -0.69 1.70 6.24 1.39 1.25 1.31 1.70 2984.14 3
## -0.69 1.70 5.67 1.95 1.25 1.31 1.70 2979.26 4
## -0.69 1.70 5.67 1.39 1.82 1.31 1.70 2976.99 5
## -0.69 1.70 5.67 1.39 1.25 1.88 1.70 2998.20 6
## -0.69 1.70 5.67 1.39 1.25 1.31 2.27 3006.84 7
## -0.53 1.87 5.83 1.55 1.41 1.47 1.14 3011.23 8
## -0.65 1.75 5.71 1.43 1.29 1.35 1.99 2984.26 9
## -0.52 1.88 5.84 1.56 1.43 0.75 1.79 2996.40 10
## -0.56 1.83 5.80 1.52 1.38 1.03 1.77 2981.02 11
## -0.53 1.86 5.83 1.54 1.41 1.19 1.44 2987.41 12
## -0.62 1.77 5.74 1.46 1.32 1.31 1.85 2977.71 13
## -0.47 1.92 5.16 1.61 1.47 1.23 1.76 2983.88 14
## -0.53 1.87 5.43 1.55 1.42 1.25 1.75 2978.32 15
## -0.59 1.80 5.49 1.49 1.35 1.57 1.70 2981.72 16
## -0.57 1.83 5.72 1.51 1.38 1.17 1.75 2976.99 17
## -0.43 1.97 5.64 0.92 1.52 1.25 1.77 2978.74 18
## -0.50 1.90 5.64 1.18 1.45 1.27 1.75 2976.32 19
## -0.58 1.81 5.94 1.22 1.36 1.31 1.73 2977.94 20
## -0.57 1.83 5.81 1.30 1.38 1.30 1.73 2976.17 21
## -0.47 1.92 5.65 1.27 1.47 1.25 1.59 2976.53 22
## -0.51 1.89 5.67 1.31 1.43 1.27 1.66 2975.55 23
## -0.35 1.32 5.72 1.32 1.59 1.24 1.73 2978.61 24
## -0.61 2.03 5.68 1.37 1.34 1.29 1.71 2975.70 25
## -0.33 1.98 5.72 1.31 0.89 1.24 1.73 2975.39 26
## -0.42 1.91 5.71 1.33 1.12 1.25 1.72 2975.05 27
## -0.41 1.88 5.67 1.14 1.26 1.40 1.68 2977.18 28
## -0.53 1.84 5.71 1.42 1.35 1.23 1.73 2975.62 29
## -0.49 1.78 5.76 1.54 1.16 1.29 1.66 2975.97 30
## -0.49 1.81 5.73 1.45 1.23 1.29 1.69 2975.32 31
## -0.40 1.86 5.57 1.46 1.19 1.26 1.67 2976.03 32
## -0.44 1.85 5.63 1.42 1.24 1.27 1.69 2975.34 33
## -0.31 1.60 5.69 1.40 1.20 1.26 1.69 2975.93 34
## -0.53 1.92 5.68 1.38 1.30 1.28 1.70 2975.26 35
## -0.39 1.81 5.65 1.34 1.18 1.34 1.66 2975.22 36
## -0.42 1.82 5.67 1.36 1.22 1.31 1.68 2975.08 37
## -0.88 1.98 5.69 1.37 1.26 1.26 1.68 2975.58 38
## -0.69 1.91 5.69 1.37 1.26 1.27 1.68 2975.25 39
## -0.54 1.81 5.69 1.45 1.03 1.30 1.73 2975.23 40
## -0.54 1.83 5.69 1.42 1.13 1.29 1.71 2975.08 41
## -0.64 1.84 5.75 1.35 1.20 1.30 1.71 2975.18 42
## -0.59 1.84 5.72 1.37 1.21 1.30 1.71 2975.07 43
## -0.62 1.88 5.65 1.30 1.20 1.29 1.72 2975.14 44
## -0.59 1.87 5.67 1.34 1.20 1.29 1.71 2975.04 45
## -0.59 1.76 5.69 1.36 1.09 1.29 1.70 2974.92 46
## -0.62 1.67 5.69 1.35 0.99 1.30 1.70 2974.99 47
## -0.41 1.72 5.69 1.36 1.09 1.31 1.73 2974.96 48
## -0.48 1.77 5.69 1.36 1.13 1.30 1.72 2974.93 49
## -0.34 1.95 5.71 1.34 1.06 1.27 1.71 2975.08 50
## -0.43 1.89 5.70 1.35 1.11 1.28 1.71 2974.98 51
## -0.61 1.86 5.72 1.36 1.07 1.26 1.74 2975.14 52
## -0.47 1.83 5.68 1.36 1.18 1.30 1.69 2974.98 53
## -0.49 1.85 5.70 1.29 1.17 1.28 1.70 2974.92 54
## -0.50 1.84 5.70 1.32 1.16 1.29 1.70 2974.92 55
## -0.40 1.83 5.66 1.32 1.08 1.28 1.71 2974.92 56
## -0.45 1.84 5.68 1.33 1.11 1.28 1.71 2974.91 57
## -0.58 1.75 5.66 1.36 1.16 1.33 1.69 2975.00 58
## -0.54 1.79 5.67 1.35 1.15 1.31 1.70 2974.94 59
## -0.40 1.77 5.71 1.36 1.07 1.30 1.70 2974.97 60
## -0.45 1.79 5.70 1.36 1.10 1.30 1.70 2974.91 61
## -0.56 1.71 5.67 1.35 1.16 1.31 1.70 2974.92 62
## -0.53 1.76 5.68 1.35 1.15 1.30 1.70 2974.90 63
## -0.54 1.75 5.69 1.33 1.08 1.29 1.71 2974.87 64
## -0.57 1.72 5.70 1.32 1.03 1.29 1.72 2974.88 65
## -0.47 1.79 5.70 1.34 1.08 1.28 1.72 2974.91 66
## -0.49 1.79 5.70 1.34 1.10 1.29 1.71 2974.89 67
## -0.53 1.81 5.69 1.32 1.09 1.28 1.70 2974.87 68
## -0.52 1.80 5.69 1.33 1.10 1.28 1.70 2974.87 69
## -0.53 1.73 5.68 1.36 1.05 1.29 1.71 2974.93 70
## -0.50 1.81 5.69 1.33 1.13 1.29 1.70 2974.89 71
## -0.41 1.83 5.69 1.32 1.12 1.29 1.71 2974.90 72
## -0.45 1.81 5.69 1.33 1.12 1.29 1.71 2974.88 73
## -0.55 1.80 5.68 1.31 1.12 1.28 1.71 2974.86 74
## -0.60 1.80 5.67 1.29 1.13 1.27 1.72 2974.87 75
## -0.58 1.75 5.70 1.33 1.11 1.29 1.70 2974.88 76
## -0.55 1.77 5.69 1.33 1.11 1.29 1.71 2974.87 77
## -0.50 1.83 5.70 1.31 1.07 1.27 1.71 2974.86 78
## -0.49 1.86 5.71 1.28 1.03 1.26 1.72 2974.89 79
## -0.55 1.81 5.69 1.31 1.10 1.28 1.70 2974.84 80
## -0.58 1.82 5.68 1.29 1.11 1.28 1.70 2974.84 81
## -0.55 1.78 5.68 1.30 1.07 1.28 1.71 2974.83 82
## -0.58 1.76 5.68 1.29 1.03 1.27 1.71 2974.82 83
## -0.64 1.77 5.69 1.29 1.06 1.28 1.71 2974.84 84
## -0.60 1.78 5.69 1.30 1.07 1.28 1.71 2974.83 85
## -0.58 1.77 5.68 1.30 1.08 1.28 1.72 2974.85 86
## -0.57 1.78 5.69 1.30 1.08 1.28 1.72 2974.84 87
## -0.57 1.81 5.68 1.28 1.05 1.27 1.72 2974.82 88
## -0.57 1.80 5.68 1.29 1.06 1.27 1.71 2974.82 89
## -0.59 1.84 5.68 1.26 1.07 1.26 1.71 2974.81 90
## -0.61 1.88 5.67 1.22 1.07 1.24 1.71 2974.82 91
## -0.59 1.80 5.69 1.27 1.02 1.27 1.71 2974.80 92
## -0.61 1.80 5.70 1.25 0.97 1.26 1.71 2974.81 93
## -0.66 1.77 5.67 1.26 1.06 1.28 1.71 2974.82 94
## -0.62 1.78 5.67 1.27 1.06 1.27 1.71 2974.81 95
## -0.59 1.77 5.68 1.28 1.00 1.26 1.72 2974.84 96
## -0.58 1.80 5.68 1.29 1.08 1.28 1.71 2974.82 97
## -0.61 1.82 5.68 1.26 1.03 1.26 1.70 2974.81 98
## -0.60 1.81 5.68 1.27 1.04 1.27 1.71 2974.81 99
## -0.59 1.82 5.67 1.25 1.03 1.26 1.71 2974.79 100
## -0.58 1.84 5.67 1.22 1.01 1.25 1.71 2974.79 101
## -0.61 1.80 5.68 1.25 1.04 1.27 1.70 2974.79 102
## -0.60 1.81 5.68 1.26 1.04 1.27 1.71 2974.79 103
## -0.62 1.86 5.68 1.23 1.06 1.26 1.70 2974.81 104
## -0.61 1.84 5.68 1.25 1.05 1.26 1.71 2974.80 105
## -0.62 1.83 5.68 1.23 1.00 1.25 1.71 2974.79 106
## -0.61 1.82 5.68 1.24 1.02 1.26 1.71 2974.79 107
## -0.62 1.79 5.68 1.24 0.99 1.26 1.71 2974.79 108
## -0.63 1.77 5.68 1.23 0.95 1.27 1.71 2974.80 109
## -0.59 1.85 5.68 1.22 0.99 1.25 1.71 2974.79 110
## -0.59 1.83 5.68 1.23 1.00 1.25 1.71 2974.79 111
## -0.60 1.83 5.68 1.21 0.99 1.25 1.71 2974.77 112
## -0.61 1.84 5.68 1.18 0.96 1.25 1.71 2974.77 113
## -0.62 1.85 5.66 1.19 1.00 1.25 1.71 2974.79 114
## -0.61 1.84 5.67 1.21 1.00 1.25 1.71 2974.78 115
## -0.61 1.82 5.67 1.20 0.95 1.25 1.71 2974.79 116
## -0.61 1.82 5.67 1.21 0.98 1.25 1.71 2974.78 117
## -0.60 1.85 5.67 1.18 0.94 1.24 1.72 2974.80 118
## -0.61 1.82 5.68 1.23 1.02 1.26 1.71 2974.78 119
## -0.59 1.82 5.67 1.21 0.99 1.26 1.71 2974.77 120
## -0.58 1.82 5.67 1.20 0.98 1.26 1.71 2974.76 121
## -0.63 1.80 5.68 1.21 0.97 1.26 1.70 2974.78 122
## -0.62 1.81 5.68 1.21 0.98 1.26 1.71 2974.77 123
## -0.59 1.86 5.67 1.18 0.99 1.25 1.71 2974.77 124
## -0.60 1.84 5.67 1.20 0.99 1.25 1.71 2974.77 125
## -0.61 1.82 5.67 1.18 0.97 1.25 1.71 2974.77 126
## -0.61 1.83 5.67 1.20 0.98 1.25 1.71 2974.77 127
## -0.60 1.84 5.67 1.17 0.94 1.25 1.71 2974.77 128
## -0.60 1.84 5.67 1.18 0.96 1.25 1.71 2974.77 129
## -0.60 1.84 5.67 1.18 0.98 1.25 1.71 2974.76 130
## -0.60 1.84 5.67 1.19 0.98 1.25 1.71 2974.77 131
## -0.59 1.83 5.68 1.18 0.95 1.26 1.71 2974.76 132
## -0.58 1.82 5.68 1.16 0.92 1.26 1.71 2974.77 133
## -0.58 1.85 5.67 1.16 0.96 1.25 1.71 2974.77 134
## -0.59 1.84 5.67 1.17 0.96 1.25 1.71 2974.77 135
## -0.59 1.83 5.67 1.19 0.98 1.26 1.71 2974.76 136
## -0.58 1.82 5.66 1.19 0.98 1.27 1.70 2974.76 137
## -0.59 1.82 5.67 1.17 0.95 1.26 1.71 2974.76 138
## -0.60 1.83 5.67 1.18 0.96 1.26 1.71 2974.76 139
## -0.59 1.82 5.67 1.18 0.97 1.26 1.71 2974.76 140
## -0.58 1.82 5.67 1.19 0.98 1.27 1.71 2974.76 141
## -0.59 1.81 5.67 1.20 0.97 1.27 1.70 2974.76 142
## -0.59 1.82 5.67 1.19 0.97 1.26 1.71 2974.76 143
## -0.57 1.83 5.67 1.19 0.97 1.27 1.71 2974.76 144
## -0.58 1.83 5.67 1.19 0.97 1.26 1.71 2974.76 145
## -0.57 1.81 5.67 1.19 0.95 1.27 1.71 2974.76 146
## -0.58 1.81 5.67 1.19 0.96 1.26 1.71 2974.76 147
## -0.59 1.82 5.67 1.17 0.95 1.26 1.71 2974.76 148
## -0.59 1.82 5.67 1.15 0.94 1.26 1.70 2974.76 149
## -0.57 1.82 5.67 1.19 0.97 1.27 1.70 2974.76 150
## -0.56 1.81 5.68 1.20 0.98 1.27 1.70 2974.76 151
## -0.57 1.81 5.67 1.20 0.99 1.27 1.70 2974.76 152
## -0.59 1.82 5.67 1.18 0.96 1.26 1.71 2974.76 153
## -0.57 1.82 5.68 1.18 0.96 1.27 1.71 2974.76 154
## -0.56 1.81 5.68 1.18 0.95 1.27 1.71 2974.76 155
## -0.56 1.82 5.68 1.18 0.96 1.27 1.71 2974.76 156
## -0.57 1.82 5.67 1.18 0.96 1.27 1.71 2974.76 157
## -0.59 1.81 5.67 1.18 0.96 1.26 1.71 2974.76 158
## -0.58 1.81 5.67 1.18 0.96 1.27 1.71 2974.76 159
## -0.58 1.82 5.67 1.17 0.96 1.27 1.70 2974.76 160
## -0.58 1.82 5.67 1.18 0.96 1.27 1.70 2974.76 161
## Assessing convergence...
## -0.47 1.82 5.68 1.18 0.96 1.27 1.71 2974.77 162
## -0.67 1.82 5.68 1.18 0.96 1.27 1.71 2974.79 163
## -0.57 1.92 5.68 1.18 0.96 1.27 1.71 2974.85 164
## -0.57 1.72 5.68 1.18 0.96 1.27 1.71 2974.78 165
## -0.57 1.82 5.78 1.18 0.96 1.27 1.71 2975.03 166
## -0.57 1.82 5.58 1.18 0.96 1.27 1.71 2975.10 167
## -0.57 1.82 5.68 1.28 0.96 1.27 1.71 2974.87 168
## -0.57 1.82 5.68 1.08 0.96 1.27 1.71 2974.83 169
## -0.57 1.82 5.68 1.18 1.06 1.27 1.71 2974.78 170
## -0.57 1.82 5.68 1.18 0.86 1.27 1.71 2974.83 171
## -0.57 1.82 5.68 1.18 0.96 1.37 1.71 2975.50 172
## -0.57 1.82 5.68 1.18 0.96 1.17 1.71 2975.34 173
## -0.57 1.82 5.68 1.18 0.96 1.27 1.81 2975.78 174
## -0.57 1.82 5.68 1.18 0.96 1.27 1.61 2975.90 175Save the model
if (fit_model) {
path <- file.path(ROOT, sprintf('notebook/models/%s', SEASON))
if (!dir.exists(path)) {
dir.create(path, recursive = T)
}
filename <- file.path(path, 'gp_la_m=200.rda')
gp_save(gp, filename)
}4 Visualise the model fit
Load the fitted model
path <- file.path(ROOT, sprintf('notebook/models/%s', SEASON))
filename <- file.path(path, 'gp_la_m=200.rda')
gp <- gp_load(filename)Make predictions in different spatial coordinates for different dates across several years.
ng <- 25
x1min <- 19.18
x1max <- 31.36
x2min <- 59.23
x2max <- 69.87
FLAG_AVERAGE_YEAR <- 0
years_all <- c(seq(2011, 2019), FLAG_AVERAGE_YEAR)
years_to_plot <- c(seq(2011, 2019, by=2), FLAG_AVERAGE_YEAR)
if (SEASON == 'spring') {
dates_to_plot <- c('1.3.', '15.3.', '1.4.', '15.4.', '1.5.', '15.5.', '30.5.')
} else {
dates_to_plot <- c('15.8.', '1.9.', '15.9.', '1.10.', '15.10.', '1.11.', '15.11.')
}
year_with_north_arrow <- FLAG_AVERAGE_YEAR
time_tol <- 8
plots_all <- list()
pr_all <- lapply(seq_along(dates_to_plot), function(i) c())
plot_index <- 1
for (year in years_all) {
average_plot <- ifelse(year == FLAG_AVERAGE_YEAR, T, F)
for (index in seq_along(dates_to_plot)) {
date_middle <- chron(paste0(dates_to_plot[index], year), format = 'd.m.y')
time <- as.numeric(as.times(date_middle))
# create (x1,x2)-grid
x1g <- seq(x1min, x1max, len=ng)
x2g <- seq(x2min, x2max, len=ng)
inputnew <- cbind(time, rep(x1g,each=ng), rep(x2g,ng) )
colnames(inputnew) <- colnames(input)
if (average_plot) {
pr <- pr_all[[index]]
signif_threshold <- 0.95 # significance level
signif <- pmax(rowMeans(pr > 0.5), rowMeans(pr < 0.5))
pr <- rowMeans(pr)
} else {
pr <- gp_draw(gp, inputnew, draws=1000)
pr_all[[index]] <- cbind(pr_all[[index]], pr)
signif_threshold <- 0.95 # significance level
signif <- pmax(rowMeans(pr > 0.5), rowMeans(pr < 0.5))
# compute the mean surface
pred <- gp_pred(gp, inputnew, transform=T)
pr <- pred$mean
# select which observations to plot
obs_ind <- input[,'time'] >= time-time_tol & input[,'time'] <= time+time_tol
obs_to_plot <- data.frame(x=input[obs_ind,'x1'], y=input[obs_ind,'x2'])
}
# create map data
mapdat <- map_data('world', regions = 'finland')
pp <- ggplot()
# where the probability is significantly different from 0.5
pp <- pp + geom_contour_filled(
data=data.frame(x=inputnew[,'x1'], y=inputnew[,'x2'], signif=signif),
aes(x=x, y=y, z=signif),
breaks=c(0.0, signif_threshold, 1.01),
alpha=0.2
)
pp <- pp + scale_fill_grey(start=1.0, end=0.5, guide='none')
# map of finland
pp <- pp + geom_polygon(data = mapdat, aes(x=long, y = lat, group = group),
fill=NA, color='black', size=0.3)
# observations
if (!average_plot)
pp <- pp + geom_point(data=obs_to_plot,
aes(x=x,y=y), color=(y[obs_ind]/n[obs_ind]>0.5)+2,
size=1, alpha=0.5, shape=16)
# prediction contours
pp <- pp + stat_contour(data=data.frame(x=inputnew[,'x1'], y=inputnew[,'x2'], pr=pr),
aes(x=x, y=y, z=pr, colour=..level..),
breaks=seq(0.1,0.9,len=7), size=0.3 )
# boundary, where probability = 0.5
pp <- pp + stat_contour(data=data.frame(x=inputnew[,'x1'], y=inputnew[,'x2'], pr=pr),
aes(x=x, y=y, z=pr, colour=..level..), breaks=0.5,
color='black', linetype=2, size=0.6 )
pp <- pp + scale_colour_gradient(limits=c(0.1,0.9), low = "red", high = "green", guide='none')
# adjust theme, xy-limits etc.
pp <- pp + theme_classic()
pp <- pp + coord_sf(crs="WGS84")
pp <- pp +
scale_x_continuous(breaks = c(20,25,30)) +
scale_y_continuous(breaks = c(60,65,70))
if ( (index == 1) && year == year_with_north_arrow ) {
pp <- pp + annotation_north_arrow(
location = "bl",
which_north = "true",
pad_x = unit(0.1, "npc"),
pad_y = unit(0.5, "npc"),
height = unit(0.2, "npc"),
width = unit(0.2, "npc"),
style = north_arrow_fancy_orienteering
)
pp <- pp + theme(
axis.line = element_blank(),
axis.title = element_blank()
)
} else {
pp <- pp + theme(
axis.line = element_blank(),
axis.ticks = element_blank(),
axis.text = element_blank(),
axis.title = element_blank()
)
}
# append
if (year %in% years_to_plot) {
plots_all[[plot_index]] <- pp
plot_index <- plot_index + 1
}
}
}
if (FLAG_AVERAGE_YEAR %in% years_to_plot) {
# change the order so that the average plot comes first
plots_year_ave <- tail(plots_all, length(dates_to_plot))
plots_year_others <- head(plots_all, length(plots_all)-length(dates_to_plot))
plots_all <- c(plots_year_ave, plots_year_others)
years_to_plot <- c(tail(years_to_plot, 1), head(years_to_plot, length(years_to_plot)-1))
}
rowtitles <- sapply(years_to_plot, function(year) ifelse(year==0, 'Average', year) )
subplot(
plots_all,
xlabs = '',
ylabs = '',
nrow = length(years_to_plot),
rowtitles = rowtitles,
coltitles = dates_to_plot,
hspace = 0.005,
vspace = 0.005,
coltitles.height = 0.15,
rowtitles.width = 0.15
)